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Abstract 

During austral summers (DJF) 1934/35, 2017/18 and 2018/19, the New Zealand (NZ) 
region (approximately 4 million km?) experienced the most intense coupled ocean- 
atmosphere heatwaves on record. Average air temperature anomalies over land were + 
1.7 to 2.1 °C while sea surface temperatures (SST) were 1.2 to 1.9 °C above average. All 
three heatwaves exhibited maximum SST anomalies west of the South Island of NZ. 
Atmospheric circulation anomalies showed a pattern of blocking centred over the Tasman 
Sea extending south-east of NZ, accompanied by strongly positive Southern Annular 
Mode conditions, and reduced trough activity over NZ. Rapid melt of seasonal snow 
occurred in all three cases. For the two most recent events, combined ice loss in the 
Southern Alps was estimated at 8.9 km? (22% of the 2017 volume). Sauvignon blanc and 
Pinot noir wine grapes had above average berry number and bunch mass in 2018 but were 
below average in 2019. Summerfruit harvest (cherries and apricots) was 14 and 2 days 
ahead of normal in 2017/18 and 2018/19 respectively. Spring wheat simulations suggested 
earlier flowering and lower grain yields compared to average, and below-average yield and 
tuber quality in potatoes crops occurred. Major species disruption occurred in marine 
ecosystems. Hindcasts indicate that the heatwaves were either atmospherically driven or 
arose from combinations of atmospheric surface warming and oceanic heat advection. 


Keywords Anthropogenic global warming - Marine heatwave - Atmospheric heatwave - 
Terrestrial ecosystems - Marine ecosystems - Crops 


1 Introduction 


Kidson (1935) described the first documented austral summer (DJF) heatwave covering the New 
Zealand (NZ) area in 1934/35, with regional temperature anomalies over land averaging + 1.7 °C 
compared to the 1981—2010 normal. At the time, this event was so unusual, almost 3 °C warmer 
than other 1930s summers; it was described as ‘remarkably warm’. Salinger et al. (2019a) 
documented the unprecedented austral summer (DJF) 2017/18 heatwave covering the NZ region. 


Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10584-020- 
02730-5) contains supplementary material, which is available to authorized users. 


AA Springer 


486 Climatic Change (2020) 162:485-506 


Regional average air (over land) and sea surface temperature anomalies were + 2.2 °C and + 1.9 °C, 
respectively. Numerous terrestrial and marine impacts persisted for the entire austral summer 
resulting in the (1) largest loss of glacier ice in the Southern Alps since 1962; (2) early Sauvignon 
blanc wine-grape maturation; and (3) major species disruption in marine ecosystems. Various 
atmospheric drivers were identified, and the event was associated with very low wind speeds, 
reduced upper ocean mixing and heat fluxes from the atmosphere to the ocean causing substantial 
warming of the stratified surface layers of the Tasman Sea. In 2018/19, New Zealand again 
experienced a very warm summer, although not as intense as the 1934/35 or 2017/18 events. 

Using the Hobday et al. (2016) definition of marine heat waves (MHW), Oliver et al. 
(2018) found a 54% increase in the number of MHW days globally since the early twentieth 
century with an increase of 3—9 days per decade in the NZ region. From two general 
circulation model (GCM) ensembles, Perkins-Kirkpatrick et al. (2018) concluded that a 
Tasman Sea MHW with the intensity of the 2017/18 event would have been virtually 
impossible without anthropogenic warming. The atmospheric blocking that was responsible 
for the prolonged period of high mean sea level pressure (MSLP) also displayed some 
anthropogenic influence, although Perkins-Kirkpatrick et al. (2018) note that this detected 
influence was less than that on the sea surface temperature. 

MHWs are caused by a range of processes at different spatial and temporal scales, from 
localized air-sea heat flux to large-scale climate drivers such as the El Niño Southern 
Oscillation (ENSO; Benthuysen et al. 2014; Heidemann and Ribbe 2019) and Southern 
Annular Mode (SAM; Thompson et al. 2011). Behrens et al. (2019) investigated mechanisms 
of MHWs in the Tasman Sea using a forced global ocean sea-ice model and Argo observa- 
tions, concluding that they are largely controlled by meridional heat transport from the 
subtropics through the interchange between the East Australian Current and the Tasman Front. 
One contributor to the increased frequency of MHWs (Oliver et al. 2018) has been regional 
warming trends. Sutton and Bowen (2019) documented a 0.1 to 0.3 °C per decade increase in 
ocean temperatures since 1981 with warming penetrating from the surface to 200 m depth 
around coastal NZ and to at least 850 m in the eastern Tasman Sea. 

This study examines the three most intense atmospheric heatwaves (AHW) and associated 
MHW for the NZ region covering the austral summers of 1934/35, 2017/18 and 2018/19. It 
diagnoses the atmospheric and oceanic drivers, impacts on marine and terrestrial ecosystems, 
including viticulture and arable cropping. Monthly to decadal atmospheric and oceanic 
mechanisms were investigated, along with an assessment of future likelihood of similar events. 


2 Methods 


Many of the methods used here were described in Salinger et al. (2019a). They are outlined 
briefly here, with new approaches described in more detail. 


2.1 Observations of atmosphere and ocean temperature 


The 22-station NZ air temperature (NZ22T) series (Salinger et al. 1992) was used to calculate 
monthly mean air temperature anomalies for 1934-2019, relative to the 1981—2010 normal. From 
daily time series, extreme statistics TX90p (percentage of days when the daily maximum temper- 
ature is above the 90th percentile), TN90p (percentage of days when the daily minimum temper- 
ature 1s above the 90th percentile) and number of summer days >25 °C averaged over New 
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Zealand during 1940-2019 were calculated as in Salinger et al. (2019a). Eight stations were 
analysed for the 1934/35 event. Monthly sea surface temperature (SST) observations were obtained 
from Extended Reconstructed Sea Surface Temperature version 5 (ERSST; Huang et al. 2017). 

For the NZ exclusive economic zone (NZEEZ), an area-weighted temperature anomaly was 
calculated from SST and land surface temperatures. The 2° x 2° ERSST product from 34° to 
48°S, and 165° to 179°E (Fig. 2g) was combined with NZ22T to produce a new temperature 
series: the NZEEZT series for 1930-2019. 

Daily SST estimates came from the NOAA /4° daily Optimum Interpolation SST analysis 
(daily OISST) on a 0.25° latitude/longitude grid spanning September 1981—July 2019 as in 
Salinger et al. (2019a). These were averaged over 160—172° E and 35-45° S (Fig. 2g). Daily 
9 am measurements of SST were obtained since 1953 at the Portobello Marine Laboratory 
(PML) in Otago Harbour, South Island, NZ. The Hobday et al. (2016, 2018) MHW definitions 
were applied to identify and characterise MHWs based on daily SST measurements from the 
daily OISST and PML datasets as in Salinger et al. (2019a). Ocean sub-surface temperature 
from GODAS (Saha et al. 2006) were averaged between 40° S and 45° S, 140° E and 150° W, 
over depths 25 to 600 m, and Argo profiles (Jayne et al. 2017) were extracted for the eastern 
Tasman Sea (160—172° E, 35—45° S). 


2.2 Atmospheric circulation 


For atmospheric circulation, monthly mean sea level pressures (MSLP) and 500-hPa geopotential 
heights were obtained from the NCEP/NCAR Reanalysis (Kistler et al. 2001) and the ERA- 
Interim reanalysis (Dee et al. 2011). Several indices were used to characterize the circulation: 
Trenberth (1976) Z1 and M1 indices and weather regimes over NZ (Kidson 2000) for the 2017/ 
18 and 2018/19 heatwaves. Z1 measures west-east (zonal) flow and M1 south-north (meridional) 
flow in the NZ region: negative Z1 is typical in blocking situations when the prevailing westerly 
to south-westerly flow is absent in the region, especially with negative M1 (northerly flow 
anomaly). For large-scale circulation and monthly to decadal modes of variability, the following 
were used: the Fogt et al. (2009) Southern Annular Mode (SAM) reconstructed index, combined 
with the Marshall (2003) SAM index, the Southern Oscillation Index (SOT) of Troup (1965) and 
for the Interdecadal Pacific Oscillation (IPO) the tripole index (Henley et al. 2015). 

A subset of past analogue (similar) three-month periods was chosen from both the ERA- 
Interim and twentieth century reanalysis (20CR, Compo et al. 2011). Analogues were selected 
for each of the three summers based on anomaly correlation and root mean-squared difference 
(RMSD) using 500 hPa anomaly fields over the NZ/Tasman Sea region. Analogues that 
exhibited anomaly correlations of at least 0.65 and RMSD of 19 geopotential metres (gpm) 
or less were selected. For 1934/35, the RMSD threshold was reduced to 16 gpm to reduce the 
number of analogue cases to a comparable level to the later summers. A total of eight 20CR 
analogues were chosen for 1934/35, while 11 and 10 ERA-Interim analogues were chosen for 
2017/18 and 2018/19 respectively. A ¢ test was used to estimate the significance of average 
SAM and SOI indices in the analogue samples chosen. 


2.3 Ocean hindcasts 
The global ocean model hindcast for this study was as described in Salinger et al. (2019a). The 


climatology for 2000—2018 was used to remove the seasonal signal. Heat content anomalies 
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were computed over the mixed layer depth, defined by a density difference of 0.01 kg/m? from 
the surface. 


2.4 Glaciers and seasonal snow 


The end of summer snowline (EOSS) time series (Chinn et al. 2012) was used to estimate 
Southern Alps mountain glacier mass balance from 1977 to 2018 for EOSSajp; (Salinger et al. 
2019a). Regression relationships were employed to calculate EOSSajp; for 1935 and 2019, 
using Hermitage Mt Cook glacier season annual mean temperature, the SAM index and 
Kidson (2000) Trough and Block regimes frequencies (Salinger et al. (2019b). EOSSprewster 
was estimated from satellite imagery and regression relations between EOSSajp, and 
EOSSpsrewster to derive a value for 2019. The methods of Chinn et al. (2012) were used to 
estimate downwasting and proglacial lake growth. 

Estimates of water stored as seasonal snow in the South Island for 2017—2018 and 2018— 
2019 were provided by the model ‘SnowSim’, available through Meridian Energy Ltd. 
(https ://www.meridianenergy.co.nz/who-we-are/our-power-stations/snow-storage/,; Garr and 
Fitzharris 1996). SnowSim calculates water stored as seasonal snow for key hydro- 
generating river catchments and is tuned to their long-term water balance. Past estimates are 
in general agreement with historical observations of snow back to 1930. Estimates for seasonal 
snow for 1934—1935 are given in Fitzharris and Garr (1995) and de Lautour (1999). 


2.5 Agriculture 
2.5.1 Horticulture 


Grapes The impact of grapevine phenology was predicted using the Grapevine Flowering 
Véraison (GFV) model (Parker et al. 2013) and harvest is defined when sugar reaches a 
concentration of 200 g/L (Parker 2012). Meteorological data (March 1947 to April 2019) were 
sourced from the Marlborough Regional Station (41.48° S; 173.95° E). Observations of yield 
component data over ten seasons (2010 to 2019) were obtained from Marlborough commercial 
vineyards (Pinot noir n = 13; Sauvignon blanc n = 34). Inflorescence numbers per metre of row 
were collected shortly after budburst in November and bunch and berry mass shortly before 
harvest. Berry number per bunch was calculated from bunch and berry mass data. 


Summer fruit Harvest dates were gathered for one variety of cherry (‘Lapins’) and two 
varieties of apricots (“Nzsummer 2’ and ‘Nzsummer3’) at the Plant and Food Research orchard 
in Clyde, Central Otago (45.20° S 169.32° E) for 2016—2019, where meteorological data were 
also obtained. 


2.5.2 Arable 


Wheat The Agricultural Production Systems sIMulator (APSIM; Holzworth et al. 2014) was 
used to assess the performance of spring wheat during the three summers 1934/35, 2017/18 
and 2018/19, against climate in the current normal period (1981—2010). Simulations were set 
for a key wheat-producing area (Lincoln, Canterbury; 43.62° S 172.47° E) by assuming that 
crops were fully irrigated and fertilized, to minimize the effects of additional yield-limiting 
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factors in the assessment. The production metrics considered were flowering time, cycle 
length, grain yield and frequency of heat stress events during the reproductive period, when 
wheat crops are sensitive to yield-damage by short periods above threshold temperatures 
(Supplementary Material). No observed crop data were available. 


Potatoes A preliminary study sought early evidence of how the 2017/18 heatwave affected 
potato production in NZ (Siano et al. 2018) at three sites: Ohakune in the central North Island 
(39.50° S, 175.45° E; 563 m above sea level (masl)) which was irrigated, Opiki, south west 
North Island (40.46° S, 175.48° E; 4 masl) which was rain-fed, and Hastings, eastern North 
Island (39.62° S, 176.73° E; 8 masl) which was irrigated. For the 2018/19 season, two locally 
bred (Ilam Hardy and Rua), and five offshore bred (Agria, Hermes, Taurus, Snowden and 
Fianna) processing cultivars were trialled at the three sites. Physiological data (net photosyn- 
thesis, transpiration rate and stomatal conductance) were measured at specific growth stages of 
the crop, while yield and tuber quality data were determined at final harvest. 


2.6 Marine ecosystems 


Impacts on marine ecosystems were evaluated from published data describing immediate 
losses of bull kelp (Durvillaea spp.) and associated community changes, and new data 
describing bull kelp community recovery 1.5 years after the 2017/18 MHW. Impacts on bull 
kelp and its community were estimated by comparing drone images and abundance surveys of 
benthic species, respectively, before and after the 2017/18 MHW from Oaro (— 42.52° S, 
173.51° E), Pile Bay (43.62° S, 172.76° E), and Moeraki (45.36° S, 170.86° E). A disturbance 
simulation experiment initiated in a bull kelp forest at Oaro in June 2017 was re-sampled in 
August 2019 to provide new data on bull kelp recovery and community change. Finally, at Pile 
Bay, we collected new before/after MHW data (where bull kelp canopy had been lost) to 
quantify changes in abundances of other canopy-forming seaweeds (see S1 and Thomsen et al. 
2019 and Thomsen and South 2019). 


3 Results 
3.1 Observations of atmosphere and oceans 
3.1.1 Surface temperatures 


The coupled ocean-atmosphere heatwaves in the NZ region during the three summer seasons 
studied here were the most intense recorded in the NZ and Tasman Sea regions in 150 years of 
land-surface air temperature records, and ~ 40 years of satellite-derived SST records (Sutton 
and Bowen 2019), as shown in Fig. la-f. 

For all three heatwaves, both land air and sea surface temperatures were 1.2° to 1.6 °C 
above the DJF 1981—2010 averages over the entire region (Salinger and Diamond 2020), from 
32° to 52°S, 150° E to 180° (Table 1 and Figs. land 2g). NZ22T anomalies (Table 1) were 
1.7 °C, 2.1 °C and 1.2 °C (Fig. la and Table 1), by far the three warmest on record (Salinger 
1979; Mullan et al. 2010). Indices of temperature extremes for NZ (Table 1 and Fig. 1b) show 
that the percentages of summer warm days above the 90th percentile were 26%, 33% and 22%, 
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Fig. 1 a New Zealand 22 station air temperatures (NZ22T) series (red smoothed) 1930-2019. b 
Extremes—TX90p, TN90p and days >25 °C averaged over New Zealand 1940-2019. c New Zealand region 
SST (34° S—48° S, 166° E-180° E) from ERSSTv5. d-f SST DJF from ERSSTv5. d. 1934/35. e 2017/18. f 
2018/19. All anomalies are from the 1981-2010 climatology period 


and percentages of warm nights above the 90th percentile were 26%, 29% and 17%. Counts of 
summer days > 25 °C averaged 22, 32 and 26 days nationwide for the three seasons. 

For the Tasman Sea and east of NZ (32°-52° S, 150°-180° E), the MHWs were 
characterised by SSTs 1.5 °C, 1.9 °C and 1.2 above average (Fig. 1d—f), the largest anomalies 
on record. All three showed a similar spatial pattern with highest anomalies to the west of the 
South Island of NZ. The largest departures from average occurred in DJF for 1934/35 and 
2017/18 but in JFM for 2018/19. 
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Table 1 Indices for the three heatwaves. NZ22T is the 22 station NZT series for surface temperature over the 
land area of New Zealand, ERSST is the ERSST version 5 for the New Zealand (NZ) Exclusive Economic Zone, 
34 to 48° S, and 165° to 179° E, and NZEEZT are NZ22Tand ERSST combined and weighted for the entire NZ 
region with anomalies in °C 


Metric 1934/35 2017/18 2018/19 
NZ22T 1.73 2.07 1.21 
ERSSTv5 1.46 1.92 n 
NZEEZ 1.48 1.92 Ty 
SAM 1.45 2.20 1.08 
SOI = 0.1 0.7 29 
IPO -0.16 = (92 0.09 
Z1 =16 = 15 =] 
M1 — 40 8 =5 
Trough 2 3 
Kidson Zonal —9 4 
Block 8 -7 
Warm days TX90p 26 33 22 
Warm nights TN90p 26 29 17 
Days = 25 °C 22 32 26 


All temperature departures are anomalies from the 1981-2010 climatology period. NZ22T mean -0.18 °C, 
standard deviation (s.d) of +0.70 °C, ERSST mean - 0.03, s.d + 0.36 °C and NZEEZT mean -0.04 s.d + 
0.37 °C. SAM the Southern Annular Mode (Fogt et al. (2009) mean 0.04, s.d 1.63, SOI the Troup (1965) 
Southern Oscillation Index mean 0.0, s.d 1.3 and the Interdecadal Pacific Oscillation (IPO) the tripolar index 
(Henley et al. 2015) mean — 0.9, s.d. 0.83, Z1 (mean — 0.4 s.d. 16) and M1 (mean — 1.7, s.d, 25) are Trenberth 
(1976) zonal and meridional indices. Kidson regimes are Trough, Zonal and Block anomalies (Kidson 2000). 
TX90p and TN90p are the percentages of days above the maximum (TX) and minimum (TN) daily 90 
percentiles (mean is 10% 1981-2010 climatology period, with days >25 °C counts >25 °C, all averaged for 
26 NZ climate stations. Mean counts days > 25 °C 15.0, s.d + 5.2 days. Italic values are significant at p < 0.05 


Applying a MHW definition (Hobday et al. 2016) to daily OISST for the eastern Tasman 
Sea (Fig. 2a, g) showed that the 2018/19 event had a similar duration but reduced intensity 
compared to the 2017/18 event. During summer 2017/18, the eastern Tasman Sea experienced 
MHW conditions for 138 days (consisting of two distinct periods of 99 days and 39 days), 
peaking as a category IV (extreme) MHW (Hobday et al. 2018) with a maximum intensity of 
4.1 °C. In comparison, the 2018/19 event lasted for 137 days, peaking as a category II (strong) 
MHW, with a maximum intensity of 2.8 °C. Nearshore surface waters at the PML followed a 
similar pattern (Fig. 2b), experiencing MHW conditions for several short (7—28 days) periods 
interspersed with cooler breaks. However, maximum anomalies during summer 2018/19 
(2.6 °C) were approximately half those observed during summer 2017/18 (5.7 °C) (Salinger 
et al. 2019a). 


3.1.2 Atmospheric circulation 


The three DJF seasons (Fig. 3a—c) show a pattern of blocking (higher than normal pressures): in 
1934/35 and 2017/18, to the east and southeast of NZ, with negative pressure anomalies 
northwest of NZ, while 2018/19 had the strongest positive pressure anomalies over the central 
Tasman Sea. The M1 and Z1 circulation indices showed northerly airflow for 1934/35, and 
easterly airflow for 1934/35 and 2017/18. Airflow was northeasterly for 2018/19. Kidson 
weather regimes showed a lack of trough types in late spring (NDJ) together with lack of zonal 
regimes and more blocking throughout the season for 2017/18. In contrast, 2018/19 had fewer 
blocking regimes, but more of the zonal regime. 
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Fig. 2 (First row) Time series of sea surface temperature (SST) climatology (1981—2011; blue), 90th percentile 
MHW threshold (orange) and summer 2017/18 to 2018/19 SSTs (black) from the a eastern Tasman Sea (160— 
172° E, 35—45° S) and b the Portobello Marine Laboratory (45.88° S, 170.5° E). The red shaded regions identify 
periods associated with MHWs from each location using the Hobday et al. (2016) definition. The red shaded 
regions identify periods associated with MHWs detected in the SST time series from each location using the 
Hobday et al. [2018] definition. (Second row) The duration of each MHW detected in the SST time series for the 
c eastern Tasman Sea and d Portobello Marine Laboratory. The red shaded region highlights MHWs detected 
between October 2017 and July 2019. (Third row) As above but showing the maximum intensity of each MHW 
detected in the SST time series for e the eastern Tasman Sea and f Portobello Marine Laboratory. (Right). g Areas 
used for SST anomalies for the New Zealand region (blue) and eastern Tamsan Sea (red) 


The 500-hPa geopotential height anomalies were consistent (Fig. 3d—f) with very 
strong blocking in the Tasman Sea extending southeast of NZ. The 1934/35 had an 
average positive height anomaly of 30 gpm west of the North Island over the north 
Tasman Sea. The 2017/18 anomaly was the most intense, 60 gpm to the south east of 
the South Island, whereas the 2018/19 anomaly was 40 gpm over the western Tasman 
Sea. The 1934/35 and 2018/19 events all exhibited ridges east and south east of the 
South Island. 

Over austral spring and summer 1934/35, the SAM was positive (Table 1). The 
SAM was also positive during 2017/18 and 2018/19. ENSO activity was weak in 
1934/35 (SOI —0.1) and in 2017/18 (SOI +0.7). The 2018/19 event was in the La 
Niña phase (SOI +2.9). The last would on average be associated with northerly 
quarter airflow anomalies in spring, and north easterly airflow anomalies over NZ 
in summer (Gordon 1986). Only in summer 2017/18 was the IPO index negative, and 
the extremely positive SAM (+ 2.2) for the season may have influenced the SSTs in 
the NZ region used to calculate the Tripole index. 


3.1.3 Analogue seasons 


The DJF 1934/35 anticyclone had a maximum anomaly of 30 gpm northwest of NZ 
(Fig. 3d). The positive anomalies covered a wide area across most of NZ and, from 
the east coast of NZ to 155° W at an average of about 10 gpm. The position 
(northwest of NZ as opposed to southeast), extent and intensity of this anticyclone 
was less than that for DJF 2017/18 season (Fig. 3e, Salinger et al. 2019a). Compared 
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Fig. 3 Atmospheric circulation patterns. DJF mean sea level pressure anomaly: a 1934/35, c 2017/18 and e 
2018/19. DJF 500-hPa Geopotential Height anomaly: b 1934/35, d 2017/18 and f 2018/19 


with the analogue cases exhibiting similar spatial patterns of circulation anomalies 
(Table 2(a)), the DJF 1934/35 SAM value (1.45) was higher than all but two of the 
analogue cases (Table 2(a)). The DJF 2017/18 season exhibited a strong SAM value 
of +2.20 (Salinger et al. 2019a) which falls in the middle of distribution of analogue 
seasons. 

The DJF 2018/19 (Fig. 3f) anticyclone anomaly reached a maximum of +40 gpm 
centred northwest of NZ. The position and extent of the positive height anomaly 
differed from 2017/18 and its intensity was considerably less. The SAM value for 
DJF 2018/19 (+ 1.08), while positive, was not statistically significant and was third 
lowest compared to the selection of analogue seasons. Therefore, the SAM likely did 
not play such a strong role in DJF 2018/19 and this may have contributed to the 
reduced warmth of the 2018/19 summer. In contrast, La Nina conditions in the 
equatorial Pacific (SOI of +7.4) would on average have been associated with warm 
SST anomalies in the Tasman Sea. For 2017/18, blocking weather types (Kidson 
2000) were most prevalent with a lack of the zonal regime and the analogue seasons 
displayed a dominance of the blocking regime. In comparison, 2018/19 zonal types 
were prevalent. In both cases, analogues exhibited a lack of the troughing regime. 


3.1.4 Ocean sub-surface temperature 


GODAS sub-surface ocean temperature patterns for DJF 2017/18 and 2018/19 for 40- 
45° S (Fig. 4) indicate very shallow positive anomalies west of the South Island, with 
a narrow band down to about 50 m east of the South Island. Positive SST anomalies 
also existed in the western Tasman Sea and into the south Pacific east of NZ. These 
were also shallow but far more intense in 2018/19 to the south east of Australia, but 
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Table 2 Detailed 500-hPa analogue results by season. These are the results of an analysis of the atmospheric 
circulation patterns were compared using anomaly correlation and root mean-square difference over the region 
135° E-140° W, 65° S—25° S (see text) compared to (a) 1934/35, (b) 2017/18 and (c) 2018/19 season. Italicized 
are significant values (p < 0.05) for SAM and SOI(/10), and 10 or 90 percentiles for Z1, M1 and Kidson regimes 


Season SAM value SOI value ZI value Mlvalue Trough  Kidson Zonal Block 
OND 1915 -0.95 —0.2 12 =26 

JFM 1935 1.00 0.4 =10 -35 

FMA 1935 — 0.50 0.3 —=23 aoa 

SON 1965 haa 1G 5 =g l 13 —14 
JFM 1966 =125 =. =o =23 = =] 10 
AMJ 1978 0.30 0.4 —40 -10 11 =] -4 
JJA 1979 2.61 =0.3 -14 — 34 — 16 6 1] 
NDJ 1981 =().50 0.5 13 = 20) 1 = 7 
AMJ 2003 2.10 -0.9 3 = 2 =] 3 5 
Mean Value 0.46 “2029 mae, sa = 0 3 
DJF 1935 1.45 -0.3 = 16 — 40 

JJA 1979 2.61 =0.3 -14 — 34 -10 3 7 
OND 1981 =().52 =(d 11 =) 15 2 11 
JJA 1985 1.06 =a = 2 = = =3 8 
DJF 1994 2.39 —0.1 l — 37 —4 = 

FMA 1999 2.70 1.2 8 =J) =13 =0 22 
MAM 1999 2.74 10 -3 —4] -18 7 J 
JAS 2005 0.02 = 8 -10 -10 9 1 
DJF 2008 2.70 1.7 =19 = = =]] 21 
DJF 2013 1.20 — 0.8 l 0 -6 5 1 
JFM 2018 1.70 4.5 = 10 0 1 E. 10 
Mean Value 1.66 0.67 -4 = 20 =li =] 10 
DJF 2018 2.20 0.7 -15 8 2 -9 8 
JFM 1982 =L 0.3 12 = 15 l me. 6 
JJA 1993 215 16 18 = 35 43 2I —8 
DJF 1995 0.77 — 0.8 -14 9 -15 LI 4 
AMJ 1999 2.32 0.7 —6 1 aoe 7 4 
MJJ 1999 1.60 0.2 0 45 = [3 16 3 
DJF 2013 1.20 — 0.8 l =22 -6 5 1 
NDJ 2015 2.52 =31 =y 3 =10 3 7 
DJF 2016 2.21 =32 =5 1 =10 5 5 
FMA 2016 1.85 3.2 1 =e 45 mae ad ma 
NDJ 2018 2.01 29 =I — 24 0 6 ~6 
Mean Value 1.49 — 0.78 =a =] -9 9 l 


DJF 2019 1.08 29 =y E 3 4 = 


not around NZ, than in 2017/18. The Argo measurements (Fig. 4c) averaged over the 
eastern Tasman Sea confirmed surface warming from December to February, peaking 
at 3 °C mean anomaly in 2017/18 and 1.5 °C in 2018/19. The anomaly was shallow, 
mainly confined to the upper 20 m when it formed, and both deepened slightly as 
they were eroded from the surface. 


3.1.5 Ocean hindcasts 


The MHW in 2015/16 which affected the region east of Tasmania (Fig. 5a) was documented by 
Oliver et al. (2017) and attributed to enhanced heat transport in the East Australian Current 
Extension (EAC-Ext). Modelled SST anomalies were intensified south of 35° S along the east 
coast of Australia and Tasmania and exceeded 1 °C. Positive SST anomalies were also present 
in the Tasman Front region, while the remaining Tasman Sea was characterised by negative 
SST anomalies. Mixed layer heat content anomalies show a pattern consistent with the SST 
anomalies, along the flow path of the EAC-Ext. where summer mixed layers are around 20 m 
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Fig. 4 Subsurface sea temperature anomalies. a, b GODAS subsurface Tasman Sea. a DJF 2017/18 and b DJF 
2018/19. c Subsurface temperature anomalies in the eastern Tasman Sea from Argo floats January 2017—April 2019 





deep. In comparison to the 2015/16 event, the 2017/18 heat wave was more intense, with SST 
anomalies above 2 °C over large parts of the Tasman Sea. The mixed layer heat content 
anomaly was positive over the entire Tasman Sea but showed a different spatial pattern 
compared to the SST anomalies, which implies differences in the driving mechanism compared 
to the 2015/16 MHW. The 2017/18 event was predominantly atmospherically driven, with low 
wind speeds reducing the vertical mixing of heat into the water column and causing a shallow 
but intense surface warming. As the surface layer warmed, the mixed layers became shallower, 
and mixed layer heat content anomalies were reduced. This differs from cases where oceanic 
heat advection is dominant and mixed layer remain shallow, as in the case of the 2015/16 and 
2018/19 MHW where SST and mixed layer heat content anomalies show similar patterns. The 
warming in 2018/19 extended from the EAC-Ext. region over the Southern Tasman Sea to the 
coastal waters of eastern NZ, and along the Chatham Rise. Each MHW event is affected by a 
different combination of both surface warming and oceanic heat advection drivers (Behrens 
et al. 2019). 

Interannual variability in Tasman Sea SST, mixed layer heat content and winds speed 
anomalies are illustrated in Fig. 5g. While the period from 2003 to 2012 was predominantly 
characterised by negative mixed layer heat content anomalies and negative SST anomalies, 
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Fig.5 Modelled SST anomalies (a—c) for November—December 2015, 2017 and 2018 in °C. The white contour 
lines show the mixed layer depths with 10 m intervals. Mixed layer heat content anomaly (d—f) for the same 
period in Joule/m2. g Timeseries of integrated or averaged anomalies over the Tasman Sea between 145° E, 175° 
E, 50° S, 30° S for November—December. Grey, red and blue bars show integrated mixed layer heat content 
anomalies in Joule, average SST anomalies in °C, and average wind speed anomalies m/s, respectively 


generally positive anomalies have been evident since then. The positive wind anomalies in 
2015, 2016 and 2017 with increased vertical mixing prevented the development of significant 
SST anomalies during spring. 


3.2 Glaciers and seasonal snow 
3.2.1 Glaciers 


Ice volume loss in the Southern Alps for the small and medium glaciers was estimated to be 
3.2 km? water equivalent (w.e.) in 1934/35, 3.6 km? in 2017/18 and 3.2 km? in 2018/19. This totals 
10 km? w.e. for the three heatwave summers, 19% of the total ice volume of the Southern Alps in 
the 1977 inventory (Chinn 2001). For the two consecutive heatwave summers, losses amounted to 
6.8 km. Total ice loss (small and medium plus 12 large glaciers) came to 8.9 km?. or an 
accumulated 22% of the 2017 volume, the largest for the entire period back to 1962 (Fig. 6a). 


3.2.2 Seasonal snow 


The 1934—1935 snow year was remarkable. Water stored as seasonal snow reached a maximum that 
was just below average at 402 mm w.e. in mid-October, based on SnowSim model estimates. Rapid 
snowmelt began in mid-November and all snow had disappeared by 11 January, the third earliest 
date for the period 1930 to 2019 (de Lautour 1999). Melt rate over this period was 6.5 mm/d w.e., 
the highest of the three summers. The earliest date for disappearance of all seasonal snow is 28 
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Fig. 6 Southern Alps ice volume and seasonal snow. a Southern Alps ice volume change (km?), between years, 
for all glaciers of the Southern Alps from 1962 to 2019. b, c Estimated water stored as seasonal snow (mm) from 
Snowsim for the period 1987-2019 


December for the 1974—75 snow year, but this was from a maximum of only 198 mm w.e., among 
the lowest since 1930. 

During the 2017—2018 snow year, the estimated water stored as seasonal snow leading up 
to August (Fig. 6b) was very low. It reached a maximum of 30% of average at 350 mm w.e. in 
late September, much earlier than usual. Rapid melt began on 18 November and from mid- 
December 2018 the snowpack was the lowest on record. By 10 January, all the seasonal snow 
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had melted, the second earliest date since 1930, with extraordinary loss of permanent glacier 
snow and ice. Melt rate over this period was 5.7 mm/d w.e. 

SnowSim estimated that maximum accumulation for the 2018—2019 snow year was close 
to average at 420 mm w.e. and occurred in late October (Fig. 6c), slightly later than normal. 
There was rapid melt from late November, but it took until 12 February for all the seasonal 
snow to disappear. Melt rate over this period was 5.0 mm/d w.e. 


3.3 Agriculture 
3.3.1 Horticulture 


Grapes Temperatures in Marlborough were above the 1981—2010 climatology for the 2017- 
2018 and 2018—2019 seasons (Table 1: supplementary Fig. S1), particularly at the key pheno- 
logical stages of inflorescence initiation (in the season before harvest), flowering and early fruit 
development (in the current season, Tables 3 and 4). Above-average temperatures at initiation and 
flowering were reflected in higher Pinot noir inflorescence number per metre of row and berry 
number per bunch (Pinot noir and Sauvignon blanc). Berry mass was reduced (Tables 3 and 4) 
supporting industry observations that the Marlborough 2019 harvest was in general less than 
anticipated (Gregan 2019). Probable environmental drivers were multiple daily maximum 
temperatures greater than 30 °C in the first 6 weeks of 2019. High temperature shock is reported 
to inhibit photosynthesis (Greer and Weston 2010) and water stress during the initial phase of 
berry development is reported to significantly reduce final berry mass (Ojeda et al. 2001). The 
GFV model simulations of flowering, véraison and harvest dates advanced since 1948 (Fig. 7) 
and the advances of last two seasons reflected the above average temperatures during spring 
(Tables 3 and 4). Despite the earlier véraison and harvest dates, mean temperatures during the 
ripening period did not increase (Fig. 1a), unlike increases observed elsewhere (Molitor and Junk 
2019). This possibly reflects the temperate climate of Marlborough and the abrupt changes in 
temperature that may occur between concurrent phenophases of vine development during the 
season (Fig. S1). 


Table 3 Key phenological stages for wine grapes, temperature and yield data. Timing of key phenological stages 
and temperatures at and between those stages 


Date Mean daily Mean daily temperature between key 
temperature at key phenology times (°C) 
phenology times (°C) 
Harvest year Flowering Véraison 20° Initiation? Flowering? 1 Sept to Flowering Véraison to 
Brix flowering to véraison 20° Brix 
Long-term 5 Dec 15 Feb 21 Mar 17.6 17.3 12.9 17.7 16.9 
average 
(1987-2017) 
2018 30 Nov 2 Feb 5Mar 17.9 18.6 13.6 19.7 19.0 
2019 2 Dec 6 Feb 12 Mar 18.7 18.5 B3 19.1 17.6 


a Initiation temperatures occur in the season before harvest, flowering temperatures are in the season of harvest. Dates 
used to estimate temperatures at this time are Dec 14 to Jan 17 and Dec 9 to Jan 9 for initiation and flowering 
respectively (Trought 2005) 
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Table 4 Sauvignon blanc and Pinot noir yield component data 

Average (2010-19) Cv? 2010-2019 = Vintage 2018 Vintage 2019 
Pinot noir 
Average bunch mass (g) 113.8 25.0 131.9 (6.5)? 76.4 (4.2) 
Average berry mass (g) 1.50 13,02 1.85 (0.04) 1.19 (0.02) 
Average berry number per bunch 39.6 13.37 46.5 (3.5) 38.9 (3.1) 
Inflorescence number per m row 24.4 19.7 25.9 (0.8) 28.1 (0.7) 
Sauvignon blanc 
Average bunch mass (g) 144.4 16.3 183.4 (3.58) 131.3 (3.63) 
Average berry mass (g) 2.05 8.15 2.37 (0.03) 1.78 (0.04) 
Average berry number per bunch 70.9 12.7 78.2 (1.3) 74.9 (2.4) 
Inflorescence number per m row 29.3 14.0 29.3 (0.6) 30.0 (0.8) 


a cv coefficient of variation 


b Numbers in brackets are the standard error of the mean 


Summerfruit Of the four seasons’ data available, September to January temperature depar- 
tures from normal during 2016-2019 were 0.0, — 0.5, +2.2 and +0.6 °C. For the cherry 
variety, 2018 and 2019 harvest dates were 13 and 3 days respectively ahead of 2016 (a normal 
season). For apricots, the two heatwave summers were 14 and 2 days ahead of normal. 


3.3.2 Arable 


Wheat APSIM-wheat simulations showed a reduction in grain yields during heatwave years 
by up to 9% compared to an estimated median historical of ~9 t/ha (Fig. 8). During heatwave 
years, there was a more frequent occurrence of shorter cycles, earlier flowering and risk of heat 
stress events throughout the reproductive phase than the historical average for Lincoln. 


Potatoes Heat and moisture stress were evident in Ohakune and Opiki (central and 
western North Island) and Hastings (eastern North Island) in 2017/18. In Opiki and 
Hastings, there were supra-optimal temperatures (>25 °C) for 54 and 60 days, 
respectively. Potato tubers from each site revealed that yield was primarily affected 
by the increase in the volume of unmarketable or defective tubers that reached up to 
85% of the total volume of tubers collected. This was due largely to the incidence of 
an array of tuber physiological defects such as enlarged lenticels, growth cracks, 
netting, malformations and second growths. 

For the 2018/19 season in Opiki and Hastings, the number of days >25 °C were 
103 and 61 days respectively (Fig. S3), with sub-optimal rainfall in Opiki (423 mm) 
(Table S2, supplementary material). In Opiki and Hastings, site average harvest index, 
total and marketable yield were reduced by up to 11.7%, 43.3% and 45.2%, respec- 
tively, with reference to the cooler environment of Ohakune. The total number of 
tubers per plant and percentage of large- and medium-sized tubers (>50 mm) de- 
clined. Dry matter content was also down by 15.7%. The elevated temperatures in 
Hastings resulted in increases in plant height and total plant leaf area, suggesting an 
enhanced dry matter partitioning to the haulm promoting vegetative growth (Levy and 
Veilleux 2007). It also increased the transpiration rate and stomatal conductance. 
Conversely, the water deficit in Opiki suppressed vegetative growth and stomatal 
conductance. These conditions led to a decrease in net photosynthesis by as much 
as 16.5%. The increase in the volume of unmarketable or defective tubers was 
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dramatic (up to 44%) which significantly reduced economic yield. The defective 
tubers exhibited physiological defects attributed to heat and moisture stress (Fig. 
S4). The most common defect was second growth in the form of heat sprouts (in- 
field sprouting), chained tubers and gemmation because of elevated soil temperatures 
and moisture stress (Hiller and Thornton 2008). Second growth was most common in 
Hastings (18.5%) and Opiki (16.8%) with extreme heat events, and lower in Ohakune 
(9.7%), which is cooler. 

The result of the trial showed location-specific adaptations (agronomic zoning) 
among the tested commercial potato cultivars. Hermes performed well in the 
drought-prone conditions of Opiki but performance was reduced in the hotter condi- 
tion of Hastings, while Snowden performed better in Hastings than in Opiki. Further 
analysis showed that ‘Taurus’, was the most stable and adaptable cultivar across test 
environments during the 2018/19 season (Fig. S5). 
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Fig. 8 Simulated physiological responses of irrigated spring wheat during three heatwave years (1934/35, 2017/ 
18 and 2018/19) in Lincoln, Canterbury, New Zealand. Dashed lines are the median (black) and average (dark- 
grey) of 30 years (1981-2010) 


3.4 Marine ecosystems 


Dramatic losses of bull kelp were reported immediately after the 2017/18 MHW, with 100% 
loss of Durvillaea poha at Pile Bay where SST exceeded 23 °C (Thomsen et al. 2019; 
Thomsen and South 2019). Cascading effects included losses of mussels and colonization of 
ephemeral seaweeds (including the invasive kelp Undaria pinnatifida). Furthermore, Undaria 
and other ephemeral seaweeds colonized plots that had lost bull kelp in Moeraki and Oaro 
(Thomsen et al. 2019; Thomsen and South 2019). New data showed that bull kelp cover was 
reduced by 60% in the undisturbed control plots in the Oaro disturbance simulation experi- 
ment. Only 1.3% of the pre-MHW juvenile bull kelp survived and no new recruits were found 
in the disturbed plots, with bull kelps being replaced by red and green turf algae. Bull kelp 
remained absent from Pile Bay and nearby reefs in Lyttelton Harbour, 1.5 years after the 2017/ 
18 MHW. After the 2017/18 MHW, plots previously inhabited by bull kelp were dominated 
by Undaria in the lower zone (43% cover) in winter months. Also, native perennial canopy- 
forming species (including Hormosira banksii, Carpophyllum maschalocarpum and 
Cystophora scalaris and Cystopohora torulosa) had recruited in 88% of the plots after the 
MHW and the loss of bull kelp (see S1 for more details). 


4 Discussion and conclusions 
Heatwaves are becoming a major impact of global warming with the Intergovernmental Panel 
on Climate Change 5th Assessment Report (Hartmann et al. 2013) indicating likely increases 


in unusually warm days and nights across most continents, and several occurrences of MHWs 
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in 2019 (Blunden and Arndt 2019). The unprecedented heatwave in the 2017/18 austral 
summer, coupled with a combined AHW/MHW event (Salinger et al. 2019a) was one 
example. Although Perkins-Kirkpatrick et al. (2018) suggests that the 2017/18 MHW would 
have been ‘virtually impossible’ without an anthropogenic influence, the 1934/35 event 
indicates a similar episode has occurred in the past which was only 0.3 °C cooler, without 
any allowance for anthropogenic global warming (AGW). Hence, it is important to examine 
similar AHW/MHWs in the NZ region in the climate record to document drivers and impacts. 

Three such austral summer events occurred—in decreasing order of magnitude 2017/18, 
1934/35 and 2018/19. The heatwaves had very similar atmospheric and oceanic footprints, 
covering all the land area, the entire central and south Tasman Sea and across to 180° E in the 
southwest Pacific. Mid-tropospheric (500-hPa) anomalies were extremely similar with 
strong blocking in the Tasman Sea extending south east of NZ. The 2017/18 event appeared 
predominantly atmospherically driven. Behrens et al. (2019) note each MHW is unique, 
either atmospherically driven, or a combination of atmospheric surface warming and oceanic 
heat advection. Trenberth et al. (2019) suggest that MHWs in the Tasman Sea region may be 
linked to heat transports from the South Pacific to the Indian Ocean north of Australia via the 
Indonesian Throughflow. Increased advection of warmer waters into the Tasman Sea is 
likely to be at the expense of a weak heat transport between the Pacific and Indian basins. As 
the inter-basin flow relaxes during El Niño years, the southward extension of the EAC is 
enhanced contributing to warming in the southern Tasman Sea. 

Projected circulation changes for the late twenty-first century (Mullan et al. 2016) show 
MSLP increases during DJF, especially to the southeast of NZ. The airflow over the country 
becomes more northeasterly, and at the same time associated with more (possibly blocking) 
anticyclones and lacking in troughs. There is also a trend towards the positive SAM resulting 
in higher MSLPs in the NZ region, but this depends on interplay with stratospheric ozone 
recovery (Arblaster et al. 2011). These are all features displayed in the 2017/18 and 2018/19 
heatwaves, with circulation regimes and their analogues exhibiting a lack of the troughing 
regime. Given that the Tasman Sea mixed layer heat content anomalies are in recent years 
have been above average, it appears likely that human-induced warming has played a 
significant role in the two recent coupled ocean-atmospheric heatwaves. 

All three heatwaves produced significant impacts on glaciers and seasonal snow and ice, 
and on terrestrial and marine ecosystems. Ice loss in small and medium glaciers has been 
estimated to range from 3.2 to 3.6 km? in each event. Across all three heatwave summers, 
there was an accumulated ice volume loss of 22% of the 2017 volume. In all three cases, 
SnowSim showed swift snowmelt commencing in mid-November and completing very 
early. Melt rates ranged between 5.0 mm/d w.e. (2018/19) and 6.5 mm/d w.e. (1934/35) 
making 1934/35 the most remarkable. 

For grapes in Marlborough, above average temperatures resulted in higher than average 
inflorescence numbers, and in 2018, berry number and bunch mass for Sauvignon blanc and 
Pinot noir wine grapes. However, 2019 berry and bunch mass were reduced, reflecting 
unusually high temperatures (over 30 °C) and vine water stress. The heat waves experienced 
in the past two growing seasons advanced the date of véraison and harvest but did not result in 
an increase in average daily temperatures during the ripening period. Harvest dates for 
Central Otago summer fruit were 2 weeks advanced in 2018 and a few days advanced in 
2019 compared with normal. 

In warm years, grain yields in wheat are reduced by the acceleration of crop development 
towards flowering and early harvest, as the crop has less time available to intercept sunlight and 
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convert it into biomass through photosynthesis. The change in flowering date also shifts the 
timing when the sensitive period to heat stress occurs, illustrating the interplay of both seasonal- 
and threshold-type damage effects in warm years (Rezaei et al. 2015). The final crop system 
response depends on various additional factors, including biotic stress and farmer’s manage- 
ment choices such as genotype selection (Teixeira et al. 2018). Nevertheless, the general 
direction of yield changes suggests greater risk to spring wheat production in heatwave years. 
For potatoes, the two heatwave years caused significant losses in the production seasons in 
terms of yield and tuber quality. 

Major species disruptions occurred in coastal marine ecosystems where bull kelp mortalities 
led to local extinctions and shifts in biodiversity. 
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